%%Solve for the BGP equilibrium of prices, wages, ouput, trade shares 

tic;

diff=10; iter=0; 



%...GIT'R DUN!...
% display('---------------------------------------------------------------');
 
while diff>10^-5;
     iter=iter+1;
w=w./w(M);

TFP=Q.^(sig-1).*Tnew;

TNT=(lamNT).*Tnew; %%Non-technology intensive productivity proportional to IP intensive one to guaranteee BGP where all countries and sectors grow at the same rate.

%% Compute prices
for i=1:M
    for j=1:M
PtempT(i,j)=(Tnew(j)*Q(j)^(sig-1))*(mbar*w(j)*distT(i,j)*(1+tau(i,j))).^(1-sig);
PtempNT(i,j)=(TNT(j)*Q(j)^(sig-1))*(mbar*w(j)*distNT(i,j)).^(1-sig);
    end
end
PT =sum(PtempT').^(1/(1-sig)); %%Price of T sector
PNT =sum(PtempNT').^(1/(1-sig)); %%Price NT sector

P=(PT./alpha).^alpha.*(PNT./(1-alpha)).^(1-alpha); %%Price index
%% Compute trade shares

for i=1:M
    for j=1:M
        pi_shareT(i,j)=Q(j)^(sig-1)*Tnew(j)*(distT(i,j)*(1+tau(i,j)))^(1-sig)*(w(j)*mbar)^(1-sig)/PT(i)^(1-sig); %%Trade share T sector
        pi_shareT_pretax(i,j)=pi_shareT(i,j)/(1+tau(i,j));
        pi_shareNT(i,j)=Q(j)^(sig-1)*TNT(j)*distNT(i,j)^(1-sig)*(w(j)*mbar)^(1-sig)/PNT(i)^(1-sig); %%Trade share NT sector
    end
end

%%Compute expenditures

YT=((pi_shareT_pretax)')^(-1)*(mbar.*w.*LT)./P'; %%Labor market clearing T sector
YNT=((pi_shareNT)')^(-1)*(mbar.*w.*LNT)./P'; %%Labor market clearing NT sector

Y=YT; %%Final production
Yw=sum(Y.*P');
%%%%%%%%%%%%%%%%%%%%%%%%%
%%With trade balance equation as the updating rule solve for wages (update wages)
%%%%%%%%%%%%%%%%%%%%%%%%%

%%First obtain royalties 

PiT= (1/(sig-1))*w.*LT; %%Profits T sector
PiNT= (1/(sig-1))*w.*LNT; %%Profits NT sector

Pi=(1/(sig-1))*w.*L; %%Total profits


%%R*beta*lam_{t+1}/lam_{t}=1 
r=(1+gg)^(1/(sig-1))/beta;  %%Discount factor (real interest rate)

for i=1:M-1
B(i)=0;
end

B(3)=-B(1)-B(2);

gp=gg/(1-sig); %%Growth prices



%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%
 %Endogenous chi
%ipr = [1 1 0.6];

    
 for i=1:M
     for j=1:M
   eps(i,j)=epsbar(i,j)*((P(i)*Ha(i,j)/Yw)^betaa);     %%Probability of adoption
 %chi(i,j) = 0.25;%(rho_chi* ((ipr(i)*Pi(i) - (P(i)*Ha(i,j))/(1+gg)))/(ipr(i)*Pi(i))) / ((1-rho_chi)*(eps(i,j)/(1-betaa))+rho_chi);       
 %chi(i,j) = (rho_chi* ((Pi(i) - (P(i)*Ha(i,j))/(1+gg)))/(Pi(i))) / (1-(1-rho_chi)*eps(i,j)*(1/r)/(1-betaa));       
 chi(i,j) = rho_chi(i);%(rho_chi(i)*((Pi(i) - (P(i)*Ha(i,j))/(1+gg)))/(Pi(i))) / (1-(1-rho_chi(i))*eps(i,j)*(1/r)/(1-betaa));       

     end
 end
 %%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%




%%Value of adopted technologies

for i=1:M
    for j=1:M
vadopter(i,j)=(1-chi(i,j))*(1-(1/r)*(1/(1+0*gp))*((1+gg)^(1/(sig-1))/(1+gg)))^(-1)*PiT(i); %%Value adopted technology to adopter
vinnovator(i,j)=(chi(i,j))*(1-(1/r)*(1/(1+0*gp))*((1+gg)^(1/(sig-1))/(1+gg)))^(-1)*PiT(i); %%Value adopted technology innovator
    end
end


 %%Guess H_{in}^a 

diffa=10; itera=0; 

     
%...GIT'R DUN!...
 %display('---------------------------------------------------------------');
 
while diffa>10^-10;
     itera=itera+1;
    
%%Probability of adoption and value of an unadopte dtechnology

for i=1:M
    for j=1:M
    %eps(i,j)=epsbar(i,j)*(1+gg)*(Ha(i,j)/P(i)/Yw)^betaa-gg;     %%Probability of adoption
   %eps(i,j)=epsbar(i,j)*((Ha(i,j)/Yw)^betaa);
    Jadopter(i,j) = (1-(1/r)*(1/(1+0*gp))*((1+gg)^(1/(sig-1))/(1+gg))*(eps(i,j)*betaa+(1-eps(i,j))))^(-1)*(eps(i,j)*(1/r)*(1/(1+0*gp))*((1+gg)^(1/(sig-1))/(1+gg))*(1-betaa))*vadopter(i,j); %%Value not yet adopted technology adopter
   Jinnovator(i,j) = (1-(1/r)*(1/(1+0*gp))*((1+gg)^(1/(sig-1))/(1+gg))*((1-eps(i,j))))^(-1)*(eps(i,j)*(1/r)*(1/(1+0*gp))*((1+gg)^(1/(sig-1))/(1+gg)))*vinnovator(i,j); %%Value not-yet adopted technology innovator (innovator chooses Hr based on this)
   %eps(i,i)=1;
    end
end

 for i=1:M
          for j=1:M
         S(i,j)=Jinnovator(i,j)^0.5*Jadopter(i,j)^0.5;     
          end
      end
      
     
%chi0=chi;
    
% for i=1:M
%     eps(i,i)=1;
% end

%%Value of innovation

for i=1:M
    for j=1:M
   vinnov(i,j)=Jinnovator(i,j)*(Tnew(j)/Tnew(i)); %%Value of an innovation (value of a patent)
    end
end
    
vinnov=sum(vinnov)';

for i=1:M
H_r(i)=P(i)^-1*(vinnov(i)*betar*lam(i)*(Yw)^(-betar))^(1/(1-betar)); %%FOC R&D spending
end


%%Get Z
ZZ=(lam'./gg).*(P'.*H_r'./(Yw)).^betar.*Tnew; %%Growth newly invente dvarieties

%%Get A_Z

for i=1:M
    for j=1:M
A_Z1(i,j)=eps(i,j)*ZZ(j)*(1+gg)/(eps(i,j)+gg); %%A(i,j)
    end
end

%A_Z2=sum(A_Z1');

for i=1:M
    for j=1:M
A_Z(i,j)=A_Z1(i,j)/Tnew(i); %%A(i,j)/T(i)
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Update H_in^a
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%Use FOC of adoption to update adoption spending

for i=1:M
    for j=1:M
LHS_Ha(i,j)=Ha(i,j)*P(i);
    end
end

RHS_Ha=gg*betaa.*(A_Z).*(1/r).*(vadopter-Jadopter).*(1/(1+0*gp)).*((1+gg)^(1/(sig-1))/(1+gg));


Za=LHS_Ha-RHS_Ha;

%% If excess < 0 country one is importing less, the increase wage in country 1 and decrease wage in 2

 % Update wages
    scl = 10; % Scale factor used to ensure that Tw > 0.
    
    THa = Ha .* (1 - scl.*Za); 
    
    % Compute distance between 'new' and 'old' wage vector.
   
     diffa = (1/scl).*max(max((abs((THa-Ha)./THa))));
   
    
    % Update wage using smoothing parameter (convex combo b/w new and old).
    smooth =  .8;
    Ha = smooth.*THa + (1-smooth).*Ha;
    
end


%%%Update wages using trade balance equation (exports T sector + exports of
%%%T sector + royalties recieves= Imports T sector + Imports NT sector +
%%%royalties paid)

  LHStemp(1:M,1) = 0;
  Lroytemp(1:M,1:M)=0;
  RHStemp(1:M,1) = 0;
  Rroytemp(1:M,1:M)=0;
  IM_T(1:M,1) = 0;
  EX_T(1:M,1) = 0;
  IM_NT(1:M,1) = 0;
  EX_NT(1:M,1) = 0;
  IBT(1:M,1) = 0;

for i=1:M
    for j=1:M
        A(i,j)=A_Z1(i,j);
    roy(i,j) =chi(i,j)*PiT(i)*A_Z(i,j); %%Royalties (fraction of profits each period) Think taht innovator sells patent to someone that then licenses the product
    end
end

for i=1:M
            for j =1:M
                if i~=j
                        IM_T(i)= IM_T(i)+ pi_shareT(i,j)*P(i)*YT(i)/(1+tau(i,j));
                         IBT(i)=IBT(i)+tau(i,j)*pi_shareT(i,j)*P(i)*YT(i)/(1+tau(i,j));
                        LHStemp(i)= LHStemp(i) +  pi_shareT(i,j)*P(i)*YT(i)/(1+tau(i,j));
                        Lroytemp(i)= Lroytemp(i) + roy(i,j);
                        LHS(i)= LHStemp(i)+Lroytemp(i);
                end
            end
        end

 for i=1:M
            for j =1:M
                if i~=j
                        EX_T(i)= EX_T(i)+ pi_shareT(j,i)*P(j)*YT(j)/(1+tau(j,i));
                        RHStemp(i)= RHStemp(i) + pi_shareT(j,i)*P(j)*YT(j)/(1+tau(j,i));
                        Rroytemp(i)= Rroytemp(i) + roy(j,i);
                        RHS(i)= RHStemp(i)+Rroytemp(i);
                end
            end
        end

%% Excess demand


Z = [(LHS'-RHS')]';

%% If excess < 0 country one is importing less, the increase wage in country 1 and decrease wage in 2

 % Update wages
    scl = 5*10^(-2); % Scale factor used to ensure that Tw > 0.
    
    Tw = w .* (1 - scl*Z'); 
    
    % Compute distance between 'new' and 'old' wage vector.
    
    diff = 1/scl*max(abs((Tw-w)./w) );
    
    % Update wage using smoothing parameter (convex combo b/w new and old).
    smooth =  0.8;
    w = smooth.*Tw + (1-smooth).*w;
    w=w./w(M);
     
    if mod(iter,100)==0
        telapsed=toc; sec=mod(telapsed,60); min=floor(telapsed/60);
        display(sprintf('\t\t SS iterations on w completed: %6.0f',iter));
        display(sprintf('\t\t\t time elapsed: %6.0f min, %2.0f sec',min,sec));
        display(sprintf('\t\t\t dif: %3.10f',diff));
    end
    
end

 %%Consumption
    
   C=Y-H_r'-sum(Ha')';